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This  paper  describes  the  development  of  a  transient  model  of  an  anode-supported,  tubular  solid-oxide 
fuel  cell  (SOFC).  Physically  based  conservation  equations  predict  the  coupled  effects  of  fuel  channel  flow, 
porous-media  transport,  heat  transfer,  thermal  chemistry,  and  electrochemistry  on  cell  performance. 
The  model  outputs  include  spatial  and  temporal  profiles  of  chemical  composition,  temperature,  veloc¬ 
ity,  and  current  density.  Mathematically  the  model  forms  a  system  of  differential-algebraic  equations 
(DAEs),  which  is  solved  computationally.  The  model  is  designed  with  process-control  applications  in 
mind,  although  it  can  certainly  be  applied  more  widely.  Although  the  physical  model  is  computation¬ 
ally  efficient,  it  is  still  too  costly  for  incorporation  directly  into  real-time  process  control.  Therefore, 
system-identification  techniques  are  used  to  develop  reduced-order,  locally  linear  models  that  can  be 
incorporated  directly  into  advanced  control  methodologies,  such  as  model  predictive  control  (MPC).  The 
paper  illustrates  the  physical  model  and  the  reduced-order  linear  state-space  model  with  examples. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

This  paper  is  the  first  part  of  a  two-paper  sequence  that  dis¬ 
cusses  the  incorporation  of  physical  models  into  advanced  process 
control.  The  present  paper  develops  transient  physical  models  for 
tubular  solid-oxide  fuel  cell  (SOFCs),  and  it  discusses  linear  model- 
reduction  strategies  near  certain  steady-state  operating  points 
(OPs).  The  companion  paper  [1  ]  concentrates  on  nonlinear  reduced 
models  and  model  predictive  control  (MPC)  design. 

The  design  and  control  of  fuel-cell  systems  can  benefit  greatly 
from  physically  based  models.  The  present  paper  is  particu¬ 
larly  concerned  with  models  that  can  be  applied  to  develop 
model-predictive  process-control  algorithms.  In  this  application, 
high-fidelity  transient  response  is  required.  The  model  considers 
flow  and  chemistry  within  a  single  anode-supported  SOFC  tube, 
as  well  as  cathode  air  that  flows  over  the  external  surfaces  of  a 
tubular  stack.  The  approach  is  based  generally  upon  earlier  works 
[2-5].  However,  to  facilitate  the  real-time  control  applications,  the 
models  must  be  computationally  fast  and  deliver  accurate  tran¬ 
sient  responses.  To  achieve  computational  efficiency,  relatively 
coarse  spatial  discretization  is  used  and  heterogeneous  reforming 
chemistry  is  modeled  with  global  reactions.  The  charge-transfer 
chemistry  is  represented  in  a  Butler-Volmer  form  that  is  derived 
from  elementary  reactions  [6]. 
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For  process-control  applications  a  pragmatic  balance  can  be 
made  between  physical  fidelity  and  computational  efficiency.  The 
primary  objective  of  the  model  is  to  assist  sensor  interpretation  and 
subsequent  actuation.  To  achieve  this  objective,  the  model  must 
incorporate  all  the  relevant  time  scales  associated  with  the  physi¬ 
cal  and  chemical  processes.  It  must  also  represent  the  relationships 
between  actuation  and  response.  However,  in  practice  the  con¬ 
trol  actuation  depends  on  feedback  from  sensors,  albeit  assisted 
by  model-based  interpretation.  In  other  words,  the  control  is  not 
based  upon  model  predictions  directly.  As  discussed  later  in  this 
paper,  complex  nonlinear  physical  models  are  further  reduced  to 
systems  of  linear  state-space  models  that  are  incorporated  into  the 
control  algorithms.  Because  the  models  that  are  incorporated  into 
the  control  algorithms  can  be  approximate,  the  underpinning  phys¬ 
ical  models  can  approximate  some  physical  attributes  in  return  for 
computational  speed. 

2.  Model  development 

Fig.  1  illustrates  a  prototype  SOFC  stack  that  is  designed  for 
a  kilowatt-scale  system.  This  stack  contains  36  anode-supported 
tubular  cells.  Fuel  enters  the  SOFC  tubes  from  below  and  exits  at 
the  top  of  the  stack.  Cathode  air  enters  radially  at  the  bottom  of 
the  stack  and  leaves  at  the  top  through  small  clearances  between 
the  tubes  and  a  top  end  plate.  The  physical  model  concentrates 
especially  on  the  flow,  chemistry,  and  thermal  behaviors  within 
the  tubes.  The  air  flow  over  the  outsides  of  the  tubes  is  mod¬ 
eled  as  a  perfectly  stirred  reactor.  The  model  accommodates  heat 
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Pf 

perimeter  between  the  fuel  channel  and  anode  sup¬ 
port  layer,  m 

Ps 

perimeter  between  the  anode  support  layer  and 
functional  layer,  m 

Pe,Pc 

perimeter  between  the  MEA  and  air  stream,  m 

P 

pressure,  Pa 

Pm2 

partial  pressure  of  H2  within  the  anode  functional 
layer,  atm 

Pf,H20 

partial  pressure  of  H20  within  the  anode  functional 
layer,  atm 

Pa,02 

partial  pressure  of  02  within  the  air  stream,  atm 

Ph2 

parameter  in  the  expression  of  i0,  atm 

Po2 

parameter  in  the  expression  of  i0c ,  atm 

Qcond 

axial  conductive  heat  flux  within  fuel  channel  flow, 
Wm-2 

Qdiff 

axial  heat  flux  within  fuel  channel  flow  from  diffu¬ 
sion,  W  m2 

^diff 

radial  heat  flux  from  the  fuel  channel  flow  into  the 
anode  from  mass  transport,  W  m-2 

^diff 

radial  heat  flux  from  the  MEA  into  the  air  stream 
from  mass  transport,  W  m-2 

Qmea 

axial  conductive  heat  flux  within  the  MEA,  W  m-2 

R 

universal  gas  constant,  J  kmol-1  K-1 

re 

outer  radius  of  anode  functional  layer,  m 

* 

outer  radius  of  fuel  channel,  m 

rP 

pore  radius,  m 

rs 

outer  radius  of  anode  support  layer,  m 

Sk 

molar  production  rate  by  surface  reactions,  kmol 
m-2  s-1 

t 

time,  s 

T 

temperature  of  the  fuel  channel  flow,  I< 

7a 

temperature  of  the  air  stream,  K 

Tm 

temperature  of  the  composite  MEA,  I< 

UAa 

overall  heat  transfer  coefficient  between  air  stream 
and  enclosure,  W  m-2  K-1 

V 

volume  of  fuel  channel  flow  control  volume,  m3 

Wk 

species  molecular  weight,  kg  kmol-1 

W 

mean  molecular  weight,  kg  kmol-1 

gas-phase  species  mole  fractions 

[Xfc] 

gas-phase  species  molar  concentrations,  kmol  m-3 

[XT] 

total  gas-phase  molar  concentrations,  kmol  m-3 

gas-phase  species  mass  fractions  within  the  fuel 
channel  flow 

gas-phase  species  mass  fractions  within  the  air 
stream  flow 

yf,k 

gas-phase  species  mass  fractions  within  the  anode 
functional  layer 

n,k 

gas-phase  species  mass  fractions  within  the  anode 
support  layer 

z 

axial  coordinate,  m 

Greek  letters 

Oi  a 

anodic  symmetric  factor  in  the  Butler-Volmer  equa¬ 
tion 

ac 

cathodic  symmetric  factor  in  the  Butler-Volmer 
equation 

<5z 

axial  length  of  fuel  channel  control  volume,  m 

A 

heat  conductivity  of  fuel  channel  flow,  W  m-1 1<-1 

Am 

heat  conductivity  of  composite  MEA,  W  m-1 1<-1 

fl 

gas-phase  viscosity,  kg  m-1  s-1 

vk 

reaction  stoichiometry  of  species  k  in  overall  charge 
transfer  reaction 

Nomenclature 

cross  sectional  area  of  fuel  channel,  m2 

As 

cross  sectional  area  of  anode  support  layer,  m2 

cross  sectional  area  of  anode  functional  layer,  m2 

^f,S 

cross  sectional  area  of  anode  support  and  functional 
layer,  m2 

as 

specific  surface  area  of  anode  support  layer,  m-1 

af 

specific  surface  area  of  anode  functional  layer,  m-1 

BS 

permeability,  m2 

Cp 

heat  capacity  of  the  fuel  channel  flow,  J  kg-1  I<-1 

Cp,m 

heat  capacity  of  the  MEA,  J  kg-1  K-1 

dp 

particle  diameter,  m 

^fc.Kn 

effective  Knudsen  diffusion  coefficients,  m2  s-1 

Dki 

binary  diffusion  coefficients,  m2  s-1 

Dh 

effective  binary  diffusion  coefficients,  m2  s-1 

Elq 

reversible  potential  between  the  anode  and  elec¬ 
trolyte,  V 

Eceq 

reversible  potential  between  the  cathode  and  elec¬ 
trolyte,  V 

Eh2 

activation  energy  for  H2  oxidation  in  the 
Butler-Volmer  equation,  J  kmol-1 

Eo2 

activation  energy  for  02  reduction  in  the 
Butler-Volmer  equation,  J  kmol-1 

EceW 

cell  potential,  V 

Erev 

reversible  cell  potential,  V 

F 

Faraday  constant,  C  kmol-1 

G° 

standard  state  Gibbs  free  energy,  J  kmol-1 

h 

enthalpy  of  the  gas-phase  mixtures,  J  kg-1 

h 

species  heat  enthalpy,  J  kg-1 

hq 

heat  transfer  coefficient  between  fuel  channel  flow 
and  MEA,  W  nr2  K-1 

hq,c 

heat  transfer  coefficient  between  air  stream  and 
MEA,  W  m-2  K-1 

*cell 

local  current  density,  A  m-2 

*0 

exchange  current  density  for  H2  oxidation,  A  m-2 

1h2 

nominal  exchange  current  density  of  H2  oxidation, 

A  m2 

*ref,H2 

parameter  i*  at  the  reference  temperature  Tref,  A 

2 

io.c 

m 

exchange  current  density  for  02  reduction,  A  m-2 

l02 

nominal  exchange  current  density  of  02  reduction, 

A  m2 

i* 

ref,02 

parameter  i*  at  the  reference  temperature  Tref,  A 

m-2 

h 

m 

gas-phase  species  mass  flux,  kg  m-2  s-1 

h 

gas-phase  species  mole  flux,  kmol  m-2  s-1 

iz,k 

axial  diffusive  mass  flux  of  k  th  species  within  the 
fuel  channel  flow,  kg  m-2  s-1 

j  r,k 

radial  mass  flux  of  k  th  species  from  the  fuel  flow 
into  the  porous  anode,  kg  m-2  s-1 

3r,k 

radial  mass  flux  of  k  th  species  from  the  MEA  into 
the  air  stream,  kg  m-2  s-1 

fr,k 

radial  mass  flux  of  k  th  species  from  the  anode  func¬ 
tional  layer  into  the  electrolyte,  kg  m-2  s-1 

3r,k 

radial  mass  flux  of  k  th  species  from  the  anode  sup¬ 
port  into  the  functional  layer,  kg  m-2  s-1 

I< 

number  of  gas-phase  species 

ne 

number  of  charge  transferred  in  the  overall  charge 
transfer  reaction 

Nu 

Nusselt  number  of  fuel  channel  flow 
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p  gas-phase  mass  density  of  fuel  channel  flow,  kg  nrr3 

pf  gas-phase  mass  density  within  anode  functional 
layer,  kg  nrr3 

ps  gas-phase  mass  density  within  anode  support  layer, 

kg  m-3 

rg  tortuosity  of  the  gas-phase 

0  g  porosity 

4>f  porosity  of  anode  functional  layer 

4>s  porosity  of  anode  support  layer 

r\  local  overpotential,  V 

Pact  local  activation  overpotential,  V 

Pact, a  local  activation  overpotential  within  the  anode,  V 

Pact.c  local  activation  overpotential  within  the  cathode,  V 

?7ohm  local  ohmic  overpotential  from  ion  conduction,  V 


cells. 


and  mass  transfer  between  the  tube  exteriors  and  the  cathode 
air. 

The  fuel-cell  model  is  written  in  terms  of  small  finite  volumes 
along  the  length  of  an  individual  tube.  Fuel  flows  within  the  tubes 
(anode  side)  and  air  circulates  outside  the  tubes  (cathode  side). 
The  fuel-flow  models  are  based  upon  a  plug-flow  approximation 
in  which  only  axial  variations  are  modeled.  That  is,  perfect  mixing 
is  assumed  across  the  radius  of  the  tube.  When  hydrocarbon  fuels 
are  used,  reforming  can  proceed  via  catalytic  chemistry  within  the 
anode  support  structure.  Charge-transfer  chemistry  proceeds  at  the 
interfaces  between  the  dense  electrolyte  and  the  composite  elec¬ 
trode  structures.  The  model  accommodates  heat  transfer  within 
the  tubular  structure  and  heat  exchange  between  the  tube  and 
gases. 

The  following  sections  present  derivations  of  the  underpin¬ 
ning  conservation  equations  for  mass  and  energy.  The  derivations 
follow  standard  practice  for  deriving  such  equations  based  upon 
appropriate  semi-differential  control  volumes.  However,  readers 


Cathode 

Dense  electrolyte 
Anode  functional  layer 
Anode  support 


Fig.  2.  Control  volume  annotated  for  the  overall  fuel-flow  continuity  equation. 


not  familiar  with  such  derivations  may  wish  to  consult  a  text  on 
chemically  reacting  flow  for  details  [7]. 


2.1.  Fuel  flow ,  overall  mass  continuity 


Fig.  2  illustrates  a  control  volume  for  a  short  tube  section.  The 
mass-conservation  equation  for  the  flowing  gases  within  tube  is 
written  as 

dp  m*-m-M  M>y 

dt  =  V  ’  (  } 

where  p  is  the  mass  density,  m*  is  the  mass  flow  rate  entering  the 
control  volume,  m  is  the  mass-flow  rate  leaving  the  control  vol¬ 
ume,  and  M  is  the  mass-flow  rate  leaving  the  control  volume  and 
entering  the  porous-anode  structure,  and  V  is  the  volume  of  the 
control  volume.  Because  the  tube  is  segmented  into  small  sections 
of  length  Sz ,  the  entering  mass  flow  rate  m*  is  equal  to  the  mass 
flow  rate  m,  which  is  leaving  the  adjacent  upstream  segment.  The 
mass  exchange  between  the  fuel  flow  and  the  anode  structure  is 
represented  as 

I< 

M  =  JjnkPfSz,  (2) 

k= 1 


where  jrk  represents  the  radial  mass  flux  of  species  kfrom  the  fuel 
flow  into  the  porous  anode  structure.  The  fuel-tube  perimeter  is 
Pf  =  277Tf  and  I<  is  the  total  number  of  species. 

A  perfect-gas  equation  of  state, 


P  1 
RT  K 

Y/k/Wk 

k=  1 


(3) 


relates  the  density,  pressure  p,  temperature  T,  and  the  mean  molec¬ 
ular  weight  W.  The  species  mass  fractions  are  Yk  and  the  molecular 
weights  are  Wk. 

By  differentiating  the  equation  of  state,  the  mass  continuity 
equation  can  be  rewritten  to  eliminate  the  mass-density  derivative 
in  favor  of  temperature  and  mass-fraction  derivatives  (pressure  is 
assumed  to  be  constant)  as 


/c-i 

-AL 


.  ■  pV  dT 

m  =  m*_M+-w+PVW^yWk 

k= 1 


J_\  9Yk 
WK)  dt  ' 


(4) 
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Fig.  3.  Control  volume  annotated  for  the  fuel-flow  species  continuity  equations. 


Fig.  4.  Control  volume  annotated  for  the  fuel-flow  thermal  energy  equation. 


The  summation  to  I<  -  1,  not  I<,  is  because  exact  species  con¬ 
servation  is  enforced  by  YK  =  1  -  Ylk=l^k-  Because  all  the  species 
mass  fractions  are  not  independent,  the  molecular  weight  WK  for 
species  I<  appears  in  the  last  term. 


2.2.  Fuel  flow,  species  continuity  equations 

Fig.  3  illustrates  species  fluxes  into  and  out  of  a  fuel  channel 
control  volume.  The  species  continuity  equation  for  the  fuel  flow 
within  the  tube  may  be  written  as 

VdJPXhl  =rhl-mk-Mk  +A,C  (; [ “k  - j{k)  .  (5) 

In  this  equation,  Yk  represents  the  mass  fractions  within  the  con¬ 
trol  volume,  which  are  the  same  as  the  mass  fractions  leaving  the 
control  volume.  The  entering  and  exiting  species  mass  flow  rates 
are 

m*k  =  m*Y*,  mk  =  mYk,  (6) 


2.3.  Fuel  flow,  thermal  energy 

As  illustrated  in  Fig.  4,  the  energy  equation  has  several  contribut¬ 
ing  terms.  In  differential  equation  form, 

^  =  (m*h*  -  mh) 

+  (“Jcond  “  ^cond)  Ac  +  “  ?diff)  Ac  ^ 

-q*mPfSz-hq(T-Tm)PfSz. 

In  this  equation  h  is  the  specific  enthalpy  and  H  =  pVh  is  the  total 
enthalpy.  The  internal  energy  time  derivative  dE/dt  is  equal  to  the 
enthalpy  time  derivative  because  the  pressure  of  the  fuel  channel 
flow  is  constant  dE/dt  =  dH/dt.  The  first  term  on  the  right-hand  side 
represents  the  energy  convected  into  and  out  of  the  control  volume. 
The  second  term  represents  thermal  conduction.  Using  Fourier’s 
law,  the  heat  flux  is  represented  as 

9cond  =  -^  (10) 

where  A.  is  the  thermal  conductivity. 

With  the  tube  discretized  into  a  series  of  control  volumes,  the 
temperature  in  each  volume  is  represented  as  T.  The  upstream  and 
downstream  conduction  are  evaluated  as 


where  Y*  are  the  entering  mass  fractions.  The  rate  of  species  mass 
exchange  between  the  anode  structure  and  the  fuel  channel  is 


=  _xT-r»  d  rd-r 

^cond  £2  ’  ^cond 


Sz 


(11) 


Mk=jr,kPfSz.  (7) 

The  axial  diffusive  mass  fluxes  across  the  upstream  and  down¬ 
stream  control  surfaces  are  represented  as  k  and  k,  respectively. 
These  axial  fluxes  are  calculated  using  “mixture  averaged”  formu¬ 
lation  [7].  The  control-surface  area  is  Afc  =  jrrf. 

By  substituting  the  overall  mass-continuity  equation  (Eq.  (1)) 
and  some  further  algebraic  manipulation,  Eq.  (5)  can  be  rewritten 
as  an  explicit  differential  equation  for  the  species  mass  fractions  as 

. K-,>- 

The  mass  fraction  of  the  I<  th  species  is  found  as  YK  =  1  - 


where  Tu  and  Td  are  the  temperatures  in  the  adjacent  upstream  and 
downstream  control  volumes.  The  third  term  in  Eq.  (9)  represents 
the  energy  that  is  carried  with  the  axial  species  diffusive  fluxes 
across  the  control  surfaces  in  the  flow  direction.  Evaluating  these 
terms  depends  upon  the  directions  of  the  species  fluxes  as 


K 


^diff  “ 


k=  1 
I< 


,  k= 1 


4%  >  0. 
ilk  <  °- 


(12) 


If  the  upstream  diffusion  flux  is  into  the  control  volume  (i.e., 
Jzk  >  0)»  t^ien  sPecies  enthalpy  is  evaluated  at  the  tempera¬ 
ture  of  the  upstream  control  volume  [i.e.,  hk{Tu)].  If  the  upstream 
diffusion  flux  is  out  of  the  control  volume  (i.e.,  j^k  <  0),  then  the 
species  enthalpy  is  evaluated  at  the  temperature  of  the  target  con¬ 
trol  volume  [i.e.,  hk{T)].  An  analogous  evaluation  pertains  to  the 
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downstream  diffusion  flux, 

I<  I< 

<4f=  E  E  03) 

fe>°  fe<° 

The  fourth  term  in  Eq.  (9)  represents  the  thermal  energy  that  is 
transferred  between  the  porous  anode  and  the  fuel  flow  as  a  result 
of  mass  exchange.  That  is, 


Qdiff  -  < 


K 


^  'jrJjhkiT ), 

k= 1 
I< 


Jr,k  >  0. 


^  'jr,l<hk(T m),  jr,k  ^  0. 
k=\ 


(14) 


Because  when  jrk  >  0  species  k  is  transported  from  the  fuel 
flow  into  the  anode  structure,  the  species  enthalpy  is  evaluated 
at  the  fuel-flow  temperature  [i.e.,  hk{T)].  When  jr  k  <  0  species  k  is 
transported  from  the  anode  structure  into  the  fuel  flow  and  species 
enthalpy  is  evaluated  at  the  membrane-electrode  assembly  (MEA) 
temperature  [i.e.,  hk{Tm)].  The  fifth  term  in  Eq.  (9)  represents  the 
thermal  energy  that  is  transferred  by  fluid  convection  between 
the  anode  surface  and  the  fuel  flow.  This  term  is  evaluated  using 
Newton’s  law  of  cooling,  where  the  heat-transfer  coefficient  frq  is 
evaluated  from  a  Nusselt-number  correlation  as 


hq 


ANu 


(15) 


For  the  laminar  flows  that  are  typical  in  small  tubular  fuel  cells, 
the  approximate  Nusselt  number  is  Nu  =  3.66.  After  substituting 
the  mass-continuity  equation  and  some  algebraic  manipulation, 
the  fuel-flow  energy  equation  can  be  written  as  a  differential  equa¬ 
tion  for  the  gas  temperature  as 


dT 

Cp~dt  ~ 


rh*  M 

7v(h*-hHWhA 

+  fend  “  qcond) 

+  feff  “  ‘Jdiff) 

k= 1 


(16) 


Because  all  species  mass  fractions  are  not  independent,  splitting 
the  time  derivative  of  specific  energy  into  temperature  and  com¬ 
position  contributions  causes  the  specific  enthalpy  hK  of  species  I< 
to  appear  in  the  last  term. 


Fig.  5.  Control  volume  annotated  for  the  bi-layer  composite  anode  structure. 


layer  and  the  dense  electrolyte  is  re.  The  average  gas-phase  compo¬ 
sitions  within  the  pore  spaces  of  the  support  and  functional  layers 
are  represented  by  mass  fractions  as  Ysk  and  Yfk,  respectively. 
Porosity  of  the  anode  support  layer  is  represented  as  0S  and  the 
specific  catalyst  surface  area  (i.e.,  effective  are  per  unit  volume)  is 
represented  as  as.  Radial  mass  fluxes  of  gas-phase  species  at  the 
interface  between  the  support  and  fuel  channel  are  represented  as 
jr  k}  radial  gas-phase  fluxes  at  the  interface  between  the  support 
and  functional  layer  are  represented  as  jsrk,  and  radial  gas-phase 
fluxes  at  the  interface  between  the  functional  layer  and  the  dense 
electrolyte  are  represented  as  k.  Axial  transport  within  the  anode 
structure  is  negligible  because  the  radial  gradients  are  significantly 
greater. 

The  species  mass-continuity  equations  within  the  anode  sup¬ 
port  layer  are  written  as 

^-t(4>sPsYs,k)  =  (as sk  +  wk,  (17) 

where  sk  are  the  molar  production  rates  of  species  via  heteroge¬ 
neous  catalytic  chemistry.  The  average  gas-phase  density  within 
the  pores  of  the  anode  support  is  ps.  The  perimeter  Pf  =  27rrf  is 
the  perimeter  between  the  fuel  channel  and  the  support  layer.  The 
perimeter  Ps  =  2tt rs  is  the  perimeter  between  the  support  layer  and 
the  functional  layer.  The  axial  cross-sectional  area  of  the  anode  sup¬ 
port  is  As  =  7t{r^  -  rf ).  The  overall  mass  continuity  equation  for  the 
support  layer,  which  is  found  by  summing  the  species  continuity 
equations  over  all  species,  is  written  as 


2.4.  Anode  bi-layer  model 

Fig.  5  illustrates  the  structure  and  notation  for  mass  balances 
within  the  porous  composite  anode  structure.  The  model  divides 
the  anode  structure  into  two  layers,  the  support  layer  and  the  func¬ 
tional  layer.  Both  layers  are  typically  the  same  composition  (e.g., 
Ni-YSZ).  The  support  layer  is  relatively  thick  (e.g.,  500-1000  pun) 
and  has  relatively  open  porosity.  There  can  be  considerable  cat¬ 
alytic  reforming  chemistry  within  the  support  layer,  but  essentially 
no  charge-transfer  chemistry.  The  functional  layer  is  relatively  thin 
(e.g.,  20-50  |jim),  with  large  three-phase-boundary  (TPB)  lengths 
and  relatively  small  particle  and  pore  spaces.  The  primary  role  of 
the  functional  layer  is  to  facilitate  charge-transfer  chemistry  near 
the  dense  electrolyte. 

The  radius  of  the  interface  between  the  support  and  functional 
layers  is  rs  and  the  radius  of  the  interface  between  the  functional 


(,8> 

k=  1 

This  continuity  equation  is  formed  assuming  no  surface  species 
are  involved  in  the  heterogeneous  reactions. 

The  species  mass-continuity  equations  within  the  anode  func¬ 
tional  layer  are  analogous  to  those  for  the  support  layer.  That  is, 

^(0fPf%)  =  (afS|£  +  J-Jlk  -  Jffk)  Wk-  (19) 

The  specific  surface  area  af  and  porosity  0f  of  the  functional 
layer  may  be  different  from  that  in  the  support  layer.  The  perimeter 
at  the  dense-electrolyte  interface  is  Pe  =  2jrre  and  the  axial  cross- 
sectional  area  is  Af  =  7r(r|  -  r^).  The  overall  continuity  equation  for 
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gases  within  the  functional  layer  is 


(20) 

k= 1 

The  heterogeneous  reforming  chemistry  is  based  upon  global 
reactions  for  steam  reforming  of  methane,  water-gas-shift  (WGS), 
and  partial  oxidation  of  methane. 


CH4  +  H2O^CO  +  3H2, 


(21) 


C0  +  H20^C02+H2,  (22) 

CH4+1o2^CO  +  2H2.  (23) 

The  rates  of  these  global  reactions  were  approximated  by  fit¬ 
ting  to  experimental  data  for  methane  reforming  as  measured  in  a 
separated  anode  experiment  [8]. 

Three  radial  mass  fluxes  must  be  evaluated  within  the  anode 
structure.  The  flux  between  the  support  and  the  fuel  flow  jrk  and 
the  flux  between  the  functional  layer  and  the  support  layer  jsr  k  are 
evaluated  using  the  dusty-gas  model  (DGM)  [9,2].  The  DGM  pro¬ 
vides  an  implicit  relationship  among  the  gas-phase  species  molar 
fluxes  through  the  porous  matrix  Jk,  molar  concentrations  [Xk], 
concentration  gradients,  and  the  pressure  p  gradient  as, 


h  [Xt]d“  d^'< 


=  -  V[Xk] 


Dk,Kn 

A 

ne  i, 


(24) 


where  [XT]  =  p/RT  is  the  total  molar  concentration,  Bg  is  the  per¬ 
meability,  and  n  is  the  mixture  viscosity.  The  mass  fluxes  jk  are 
related  simply  to  the  molar  fluxes  Jk  as  jk  =  WjJk.  Knudsen  dif¬ 
fusion  represents  mass  transport  assisted  by  gas-wall  collisions. 
The  Knudsen  diffusion  coefficients  depend  upon  the  porous-media 
microstructure,  including  porosity,  average  pore  radius  rp,  and  tor¬ 
tuosity  Tg.  The  effective  binary  and  Knudsen  diffusion  coefficients 
Dkl  and  Dk  Kn  can  be  evaluated  as 


=  DLcn  =  i 


2  rP0g 

3  r E 


8  RT 
*Wk' 


(25) 


The  binary  diffusion  coefficients  Dki  and  the  mixture  viscosities 
fi  are  determined  from  kinetic  theory  [7].  The  permeability  can  be 
evaluated  from  the  Kozeny-Carman  relationship  as 


Bg  = 


72rg(l  -  4>g)2  ’ 


(26) 


where  dp  is  the  particle  diameter.  Further  details  of  the  DGM  and  its 
numerical  implementation  can  be  found  in  Zhu,  et  al.  [2].  The  two 
fluxes  jr  k  and  jsr  k  are  evaluated  based  on  the  transport  properties 
of  the  anode  support. 

The  mass  fluxes  at  the  interface  between  the  anode  functional 
layer  and  the  dense  electrolyte  k  are  determined  by  the  charge- 
transfer  chemistry 


Jr,k  = 


^lAcell 

neF 


Wk, 


(27) 


where  icen  is  the  current  density  and  F  is  Faraday’s  constant.  The 
variables  vk  and  ne  are  the  stoichiometric  coefficient  for  species  k 
and  number  of  electrons  transferred  in  the  overall  charge  transfer 
reaction. 


2.5.  MEA  energy  balance 


Fig.  6.  Control  volume  annotated  for  the  thermal  balance  in  the  MEA,  which  is  the 
entire  tube  wall. 


Fig.  6,  heat  enters  and  exits  the  MEA  through  conduction,  mass 
transport,  and  convection.  The  energy  balance  for  the  MEA  is  writ¬ 
ten  as 


( 1  —  0m  ) 


^(PmCp^Tm) 


-  (au  -ad  )  — 

—  ^mea  4mea )  ^ 


dt 

I< 

Pf  I I  Pe  -r  u 

+/^sY^Jr-khk~  ArJr<°2h°2 


k= 1 


^f,s 


T^~hq(T  —  Tm)  —  hq)C(Tm  —  Ta. 
—  ^celdcelb 


(28) 


where  Afs  is  the  combined  cross  sectional  area  of  the  anode  sup¬ 
port  and  functional  layer.  The  first  term  on  the  right-hand  side 
represents  axial  thermal  conduction  through  the  solid  materials 
that  comprise  the  MEA.  The  upstream  and  downstream  conduction 
terms  are  evaluated  as 


9mea  — 


T  Tu 
_  m  -  jrn 

8z 


9mea 


—  — 


yd  t 
1m~  1  m 

Sz 


(29) 


where  and  are  the  MEA  temperatures  in  the  adjacent 
upstream  and  downstream  control  volumes  and  Am  is  the  effec¬ 
tive  thermal  conductivity  of  the  MEA.  The  MEA  porosity  0m  is  an 
area  average  of  anode-support  and  anode-functional  layers  poros¬ 
ity.  The  second  term  on  the  right-hand  side  of  Eq.  (28)  represents 
the  energy  that  is  transported  via  species  flux  from  the  fuel  stream 
to  the  anode  support  layer.  The  third  term  represents  the  energy 
that  is  transported  via  species  flux  from  the  cathode  to  the  exte¬ 
rior.  As  discussed  in  the  context  of  the  fuel-flow  energy  balance 
(e.g.,  Eq.  (14)),  the  temperature  at  which  the  enthalpy  is  evaluated 
depends  upon  the  direction  of  the  species  flux.  In  a  similar  way, 
energy  exchange  between  the  cathode  and  the  exterior  is  evaluated 
as 


'  I< 

^  Jr,k^kO m)> 


Qdiff  -  < 


k= 1 
I< 


<  k= 1 


jc.  >0 

Jr,k 

f,  <o 

Jr,k 


(30) 


The  entire  thickness  of  the  MEA  is  represented  by  a  single  tem¬ 
perature  Tm,  which  varies  temporally  and  axially.  As  illustrated  in 


Because  when  jcr  k  >  0  species  k  is  transported  from  the  cathode 
into  the  surrounding  air  stream,  the  species  enthalpy  is  evaluated 
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at  the  MEA  temperature  [i.e.,  hk(Tm)].  When  j^k  <  0  species  k  is 
transported  from  the  surrounding  air  stream  and  species  enthalpy 
is  evaluated  at  the  air  stream  temperature  [i.e.,  hk{Ta)]. 

The  fourth  and  fifth  terms  in  Eq.  (28)  represent  convective  heat 
transfer  between  the  MEA  and  adjacent  gas  flows.  The  heat-transfer 
coefficients  on  the  exterior  cathode  side  hqtC  is  an  empirical  param¬ 
eter.  The  final  term  in  Eq.  (28)  represents  the  electrical  energy 
produced  by  the  cell,  which  does  not  contribute  to  the  thermal 
energy  balance.  Assuming  that  the  dense  electrolyte  and  cathode 
are  very  thin,  the  cathode  perimeter  is  Pc  =  Pe  =  2jzre. 

2.6.  Cathode  air  conservation  equations 


current  density,  Eq.  (34)  can  be  used  to  determine  the  local  current 
density  for  each  discrete  control  volume  along  the  length  of  the 
SOFC  tube.  Because  the  model  determines  the  composition  in  the 
anode  functional  layer,  there  is  no  need  to  include  concentration 
overpotentials. 

The  reversible  voltage  is  calculated  based  on  the  global  charge- 
transfer  reactions  at  the  anode  and  cathode  side  being  in  local 
equilibrium.  On  the  cathode,  oxygen  from  the  air  stream  is  reduced 
by  reaction  with  electrons,  producing  an  oxygen  ion  in  the  elec¬ 
trolyte  as 

lo2(g)  +  2e-(c)^02-(el).  (35) 


The  model  is  based  upon  an  assumption  that  the  cathode  air  that 
circulates  around  the  outsides  of  the  tube  stack  behaves  as  a  per¬ 
fectly  stirred  reactor.  That  is,  the  air  temperature  and  composition 
may  vary  temporally,  but  are  uniform  spatially.  The  rate  at  which 
oxygen  is  transferred  from  the  circulating  air  to  the  SOFC  cath¬ 
odes  depends  upon  the  cell  operating  conditions.  Heat  is  exchanged 
between  the  tubes  and  the  air  by  convection  and  by  the  energy 
associated  with  the  oxygen  mass  transfer. 

The  species  continuity  equation  for  the  cathode  air  is  written  as 

Va^^=rh*y*fc-rhaya,k,  k  +  02,  (31) 

where  m*  is  the  mass  flow  rate  of  the  inlet  air  and  Y*  are  the  inlet 

a  a  ,k 

species  mass  fractions  (usually  only  02  and  N2 ).  The  mixture  within 
the  volume  Va  of  cathode  space  and  the  exhaust  flow  are  assumed  to 
have  the  same  composition,  Yak.  The  oxygen  mass  fraction  is  deter¬ 
mined  from  Ya  q2  =  1  -  Ylk=j  o2^a-fc‘  overa^  mass-continuity 
equation  is 

(32) 


where  Ma  is  the  oxygen  mass-consumption  rate  by  the  SOFC  stack. 
The  cathode-air  energy  equation  is 


,,  d(pe) 
Va  dt 


m*h*  -  maha  -  UAa(Ta  -  TsheU) 


+Pe 

+^e 


hq,c(Tm  —  Ta)dz 

Jr,  02ho2dz. 


(33) 


The  first  two  terms  on  the  right-hand  side  represent  the  energy 
associated  with  the  air  flow  into  and  out  of  the  cathode-air  com¬ 
partment.  The  third  term,  which  represents  heat  transfer  from  the 
cathode  air  to  the  surrounding  enclosure,  uses  an  overall  heat  trans¬ 
fer  coefficient  UAa.  The  fourth  term  is  an  integral  that  represents  the 
convective  heat  transfer  along  the  length  L  of  an  SOFC  tube  to  the 
surrounding  air.  The  final  term  represents  the  energy  associated 
with  oxygen  mass  exchange  between  the  SOFC  cathode  and  the 
surrounding  air. 


2.7.  Electrochemistry 

The  first  step  in  modeling  the  electrochemistry  is  to  determine 
the  cell  voltage  £cen.  The  cell  voltage  is  equal  to  the  reversible  volt¬ 
age  minus  the  sum  of  various  overpotentials 

ECell  —  ^rev  —  Pact, a  +  Pact.c  ~  Pohm>  (34) 

where  Erev  is  the  reversible  voltage  (cell  voltage  when  no  net  cur¬ 
rent  is  produced).  The  ohmic  overpotential  pohm  is  mainly  from  ion 
conduction  through  the  electrolyte.  There  is  an  activation  overpo¬ 
tential  to  drive  charge-transfer  chemistry  at  the  cathode  pact,c  and 
the  anode  pact,a-  Because  each  of  the  overpotentials  depends  upon 


The  oxygen  ion  is  transported  across  the  dense  electrolyte 
toward  the  anode  functional  layer.  Within  the  anode  functional 
layer,  hydrogen  from  the  fuel  stream  reacts  with  the  oxygen  ion 
and  delivers  electrons  into  the  anode  phase  as 

H2(g)  +  02-(el)  -  H20(g)  +  2e-(a).  (36) 


The  overall  reversible  potential  can  be  calculated  by  setting  the 
electrochemical  energy  change  for  reactions  (35)  and  (36)  to  zero. 
The  result  is 


AG°  RTm 

Erev  2F  +  "2F  n 


Vh -AS, 

Pf,H20 


(37) 
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f,H20 
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where  pk  f  and  pk  a  are  the  partial  pressures  (evaluated  in  atm)  of 
species  k  in  the  anode  functional  layer  and  air  stream,  respectively. 
The  standard  state  Gibbs  free  energy  change  G°  is  only  a  function 
of  temperature  and  not  concentration.  Thus,  the  first  term  in  Eq. 
(37)  accounts  for  temperature  effects.  The  last  term  accounts  for 
composition  and  pressure  effects. 

Butler-Volmer  equations  are  used  to  determine  the  local  cur¬ 
rent  density  i  at  the  anode  and  cathode  side  for  a  given  activation 
overpotential 


The  exchange  current  density  is  a  function  of  the  local  con¬ 
centrations  and  temperature  at  the  TPBs.  The  global  anodic  and 
cathodic  symmetry  factors  aa  and  ac  do  not  necessarily  sum  to 
one.  For  a  given  charge-transfer  reaction  mechanism,  the  functional 
form  of  the  exchange  current  density  and  global  symmetry  factors 
are  derived  by  assuming  one  elementary  reaction  is  rate  limiting 
and  that  all  other  reactions  are  in  partial  equilibrium. 

The  assumed  hydrogen  oxidation  mechanism  and  rate  limiting 
step  are  taken  from  Zhu  et  al.  [2].  By  assuming  the  symmetry  fac¬ 
tors  for  the  elementary  rate  limiting  step  are  both  aa=ac  =  0.5, 
the  functional  form  of  the  exchange  current  density  for  the  elec¬ 
trochemical  oxidation  of  hydrogen  is 


.  (Pf3G/Pii2)1/4(Pf.H2o)3/4 

10  'H2  l+(Pf,H2/Pf,J1/2 


(40) 


where  i^2  and  p^  are  only  functions  of  temperature.  Also,  the 
global  symmetry  factors  for  hydrogen  oxidation  are  aa  =  1.5  and 
ac  =  0.5.  For  the  cathode,  the  mechanism  for  the  electrochemical 
reduction  of  oxygen  is  also  taken  from  Zhu  et  al.  [2].  Based  on  the 
assumed  elementary  rate  limiting  step  having  symmetry  factors 
equal  to  0.5,  the  exchange  current  density  for  the  cathode  side  is 


'0,c  -  'o2 


(pa,o2/Po2)1/4 
1  +(Pa,o2/P52)1/Z 


(41) 
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The  global  symmetry  factors  for  the  electrochemical  reduction 
of  oxygen  are  aa  =  0.5  and  ac  =  0.5.  The  nominal  exchange  current 
densities  i*  depend  on  temperature  in  the  following  way 


it  =  i’ref  e  exp 


(  Ee 

1  ^ 

\  R 

(t 

TrefJJ 

(42) 


where  “e”  represents  the  electrochemical  oxidation  of  oxygen  or 
hydrogen.  The  parameter  i*efe  is  the  nominal  exchange  current 
density  at  the  reference  temperature  Tref.  It  is  important  to  note 
that  all  partial  pressures  in  Eqs.  (40)  and  (41 )  must  be  evaluated  in 
atmospheres. 


3.  Implementation 

Spatial  derivatives  are  evaluated  using  the  finite  volume 
method.  The  resulting  governing  equations  from  a  set  of 
differential-algebraic  equations  (DAEs)  [10].  The  model  is  imple¬ 
mented  in  C++  and  linked  to  Sundials[11]  and  Cantera[12].  The 
discretized  equations  are  solved  using  the  DAE  solver  IDA  con¬ 
tained  within  Sundials.  Based  on  mesh  studies,  an  axial  mesh  of 
100  points  accurately  captures  spatial  variation  along  the  tube. 
For  each  axial  mesh  point,  there  is  an  anode  support  layer,  anode 
functional  layer,  and  fuel  channel  control  volume.  The  chemically 
reacting  flow  software  Cantera  is  utilized  to  calculate  fluxes  within 
the  anode  according  to  the  DGM.  Cantera  is  also  used  to  evaluate 
thermodynamic  properties,  transport  properties,  and  net  rates  of 
production  from  thermal  chemistry. 


4.  Example  model  results 

A  specific  anode-supported  SOFC  tube  is  used  to  illustrate  the 
model.  The  15  cm  long  tube  with  an  outside  diameter  of  1  cm., 
is  typical  of  a  tube  that  may  be  used  in  a  sub-kilowatt  tubular 
stack.  Table  1  lists  other  model  parameters,  describing  physical 
and  electrochemical  characteristics.  Most  of  these  parameters  are 


Table  1 

Model  parameters  for  a  single  tube.  Many  of  the  parameters  are  from  Zhu  et  al.  [4]. 


Parameter 

Value 

Thermal  properties 

MEA  thermal  conductivity,  A.m  (W  irr 1  K-1 ) 

10.5 

MEA  porosity,  </>m 

0.35 

MEA  heat  capacity,  cPim  (J  kg-1  K_1 ) 

533 

MEA  density,  pm  (kg  rrr3) 

7000 

Nominal  heat  transfer  coefficient,  hq,c  ( W  m-2  I<-1 ) 

120 

Nominal  overall  heat  transfer  coefficient,  UAa  (W  I<-1 ) 

0.85 

Nusselt  number,  NuD  =  hqPf/(nk) 

3.66 

Anode  support  layer  properties 

Thickness,  Ls  (/im) 

950 

Porosity,  0S 

0.35 

Tortuosity,  r 

3.5 

Pore  radius,  rp  (/xm) 

0.5 

Mean  particle  diameter,  dp  (/im) 

2.5 

Specific  catalyst  area,  as  (m-1 ) 

1.08e5 

Anode  functional  layer  properties 

Thickness,  Lf  (/x  m) 

50 

Porosity,  0f 

0.35 

Specific  catalyst  area,  af  (m_1 ) 

2.08e5 

Nominal  exchange  current  density,  i*2  (A  m-2) 

2.07 

Activation  energy,  Eh2  (kj  mol-1 ) 

120 

Reference  temperature,  Tre f  (°C) 

800 

Electrolyte  ohmic  resistance  R  =  R0T  exp (Eion/(RT)) 

Activation  energy,  Eion  (kj  mol-1 ) 

80 

Resistance  prefactor,  R0  (C2m2  I<-1 ) 

5.55e-13 

Cathode  properties 

Nominal  exchange  current  density,  i*2  (A  m-2) 

0.68 

Activation  energy,  Eo2  (kj  mol-1 ) 

130 

Reference  temperature,  Tre f  (°C) 

800 

taken  from  Zhu  et  al.  [4].  The  effective  thermal  conductivity  of  the 
composite  MEA  (i.e.,  tube  wall)  is  approximated  as  Am  =  10.5  W 
m-1  K-1  [13].  The  nominal  exchange  current  densities  =  2.07  A 
cm-2  and  i*2  =  0.68  A  cm-2  were  adjusted  such  that  model  pre¬ 
dicted  power  densities  match  those  of  typical  small  scale  systems 
(around  0.3  W  cm-2  for  fuel  utilization  of  around  85%). 

For  purposes  of  illustrating  the  model,  the  inlet  fuel  composi¬ 
tion  is  38%  H2, 3%  H20, 1%  CH4, 19%  CO,  0.3%  C02,  and  38%  N2.  This 
is  the  equilibrium  mixture  for  methane  and  air  at  800°C  with  a  sto- 
chiometric  ratio  to  partial  oxidation  (i.e.,  oxygen  to  carbon  ratio  of 
0.5).  A  small  amount  of  steam  (3%)  has  been  added  to  the  equi- 
lbrium  mixture.  Fuel  enters  the  tube  at  800  °C  and  atmospheric 
pressure.  The  cathode  air,  which  flows  around  the  outsides  of  the 
tubes,  enters  the  stack  at  550°C.  Both  ends  of  the  tubes  conduct 
heat  to  manifolds  at  a  fixed  temperature  of  800°C. 

4.1.  Steady -state  results 

Fig.  7  illustrates  steady-state  spatial  profiles  for  composition, 
current  density,  temperature  and  fuel  velocity  with  the  cell  oper¬ 
ating  at  £cell  =  0.72  V.  The  fuel  inlet  velocity  is  38  cm  s-1,  and  the  air 
flow  rate  (per  tube)  is  33.5  mg  s-1 .  The  total  current  and  power  pro¬ 
duced  by  a  single  tube  are  23.1  A  and  16.6  W,  respectively,  which 
corresponds  to  current  and  power  densities  of  0.49  A  cm-2  and  0.35 
W  cm-2.  Under  these  conditions,  the  cell  achieves  an  efficiency  of 
49.9%  and  fuel  utilization  of  89.9%. 

The  fuel  composition  profiles  (Fig.  7)  are  primarily  driven  by 
the  electrochemical  consumption  of  hydrogen.  The  ratio  of  hydro¬ 
gen  to  steam  decreases  along  the  length  of  the  cell  due  to  current 
production.  The  heterogeneous  WGS  process  (Eq.  (22))  remains 
near  equilibrium.  Thus,  carbon  monoxide  and  steam  react  rapidly, 
replenishing  the  electrochemically  consumed  hydrogen.  Overall, 
these  processes  consume  hydrogen  and  carbon  monoxide,  produc¬ 
ing  steam  and  carbon  dioxide.  The  small  amount  of  methane  (1%) 
entering  the  tube  is  quickly  reformed  by  steam  (Eq.  (21))  to  pro¬ 
duce  CO  and  H2.  Because  this  reaction  increases  moles,  the  nitrogen 
mole  fraction  decreases  slightly.  The  methane  steam  reforming  also 
produces  a  slight  maximum  in  the  velocity.  The  peak  in  the  tem¬ 
perature  profile  is  the  result  of  internal  heat  generation  and  fixing 
the  end  temperatures. 

Fig.  7  illustrates  some  interesting  effects  of  temperature  and  fuel 
composition  on  local  current  density.  Near  the  tube  inlet,  the  cur- 
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Fig.  7.  Steady-state  profiles  for  temperature,  fuel  channel  composition,  current  den¬ 
sity,  and  fuel  velocity  along  the  length  of  the  cell.  The  cell  voltage,  fuel  inlet  velocity, 
and  air  mass  flow  rate  are  0.72  V,  38  cm  s-1  and  33.5  mg  s-1,  respectively. 
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rent  density  increases.  The  current  density  reaches  a  maximum  of 
0.74  A  cm-2  about  2  cm  downstream  of  the  inlet.  Following  the 
maximum,  the  current  density  decreases  monotonically  to  around 
0.17  A  cm-2  at  the  fuel  outlet.  The  reversible  potential  decreases 
along  the  tube  as  hydrogen  is  electrochemically  consumed  and 
steam  is  produced,  tending  to  reduce  the  local  current  density. 
However,  the  current  density  also  depends  upon  the  local  MEA 
temperature,  which  is  maximum  at  around  3.5  cm  from  the  tube 
entrance.  These  competing  effects  lead  to  the  maximum  current 
density  occurring  between  the  fuel  inlet  and  maximum  MEA  tem¬ 
perature. 

4.2.  Transient  simulations 

The  model  is  designed  to  run  simulations  with  time-varying 
inputs.  For  instance,  the  cell  voltage,  fuel  inlet  velocity,  and 
air  flow  rate  can  be  adjusted  to  meet  varying  power  demands. 
Consequently,  the  spatial  profiles  (composition,  velocity,  cur¬ 
rent  density,  etc.)  are  transient.  However,  the  sensors  are 
assumed  to  be  positioned  to  measure  tube-exhaust  conditions. 
Thus,  the  controller  is  more  concerned  with  the  temporal 
variation  of  the  outlet  flows,  average  MEA  temperature,  and 
total  current  production  than  the  spatial  variations  within  the 
tube. 

Consider  the  transient  behavior  that  results  in  going  from  a  low 
current  demand  to  a  high  current  demand.  One  second  into  the  sim¬ 
ulation,  the  current  demand  goes  from  1 4.5  A  to  30  A.  At  the  start  of 
the  simulation,  the  cell  is  operating  in  steady  state  with  a  cell  volt¬ 
age  of  Eceii  =  0.78  V,  producing  1 1.3  W.  The  inlet  fuel  velocity  is  25 
cm  s-1  and  air  flow  rate  is  33.5  mg  s-1 .  Under  these  conditions,  the 
cell  is  operating  at  52.3%  efficiency  and  a  fuel  utilization  of  87.2%. 
Figs.  8  and  9  illustrate  transient  response  on  short  (seconds)  and 
long  (minutes)  time  scales,  respectively. 

One  second  after  the  start  of  the  simulation,  the  current  demand 
is  increased  from  14.5  A  to  30  A.  In  an  attempt  to  achieve  the  new 
current  demand,  the  cell  voltage  is  reduced  from  Ecell  =  0.78  V  to 
ECeii  =  0.69  V.  Lowering  the  cell  voltage  further  could  permanently 
degrade  the  tube.  The  voltage  drop  causes  a  sharp  current  increase 
from  14.5  A  to  around  25  A.  The  electrochemical  charge-transfer 
rate  responds  instantly  to  the  change  in  cell  voltage.  However,  the 
suddenly  increased  current  production  increases  fuel  consumption, 
causing  a  subsequent  decrease  in  current  (Fig.  8).  Thus,  the  new  cur¬ 
rent  demand  of  30  A  is  not  met  with  only  a  change  in  cell  voltage. 
Also,  the  fuel  utilization  approaches  100%  (i.e.,  very  low  H2  in  the 
exhaust),  which  is  an  undesirable  condition.  To  meet  the  current 
demand  of  30  A  and  offset  fuel  depletion,  the  fuel  inlet  velocity  is 
increased  from  24.7  cm  s-1  to  51 .4  cm  s-1  over  a  half-second  ramp 
beginning  at  1.2  s.  The  increased  fuel  flow  causes  the  current  to 
increase  to  around  28  A  and  the  H2  in  the  exhaust  increases  to 
around  7%  by  around  2  s.  Based  upon  the  characteristic  flow  and 
diffusion  times,  the  response  time  for  rapid  changes  in  operating 
conditions  is  around  1  s. 

Beginning  at  2  s  after  the  start  of  simulation,  the  current 
increases  on  a  shallow  linear  ramp  toward  and  above  the  desired 
value  of  30  A.  With  the  voltage  and  inlet  velocity  held  fixed,  the 
fuel  consumption  increases  leading  to  a  decrease  in  the  H2  exhaust 
mole  fraction  (Fig.  8).  During  the  first  10  s,  the  tube  temperature, 
fuel-exhaust  temperature,  and  cathode  air  temperature  all  increase 
slightly.  The  temperature  rise  is  the  result  of  increasing  polarization 
losses  as  the  cell  power  output  increases. 

Fig.  9  shows  cell  behavior  on  a  longer  time  scale  extending  to 
160  s.  After  60  s,  the  current  continues  to  rise  above  the  demanded 
value  of  30  A  due  to  a  continued  rise  in  temperatures.  The  increased 
current  and  increased  temperature  cause  increased  fuel  consump¬ 
tion.  If  the  tube  temperature  becomes  sufficiently  high  for  a  long 
period  of  time,  then  the  tube  might  be  permanently  damaged.  To 


Fig.  8.  Total  current,  exhaust  H2  mole  fraction,  and  temperatures  as  functions  of 
time,  following  a  transient  from  a  relatively  low  power  to  high  power  demand.  The 
change  is  effected  by  reducing  the  cell  voltage  from  0.78  V  to  0.69  V  one  second  after 
the  start  of  simulation.  The  velocity  of  the  fuel  inlet  is  increased  from  24.7  cm  s-1  to 
51.4cms_1  over  a  half-second  ramp  beginning  at  1.2  s  after  the  start  of  simulation. 
The  responses  on  the  short  time  scale  are  characterized  by  the  fluid-flow  dynamics 
that  are  on  the  order  of  seconds. 


stabilize  the  current  production  at  30  A  and  lower  the  average  tube 
temperature,  at  around  80  s  the  cathode  air  flow  rate  is  increased. 
The  cooling  effect  of  the  cathode  air  serves  to  decrease  tube  and 
flow  temperatures,  although  on  relatively  long  time  scales.  The 
average  tube  temperature  decreases  and  stabilizes  at  830°  C.  The 
tube  current  stabilizes  at  29.9  A.  As  the  temperatures  decrease,  the 
fuel-consumption  rate  decreases  (H2  in  exhaust  increases)  and  the 
total  current  decreases  slightly.  At  the  end  of  simulation  (160  s),  the 
cell  is  operating  with  an  efficiency  of  45.9%  with  a  fuel  utilization 
of  85.9%. 

This  simulation  illustrates  the  need  for  control  of  fuel  cell  sys¬ 
tems  to  quickly  meet  varying  current  demands  and  not  violate 
constraints.  The  desired  higher  current  demand  was  only  met  and 
stabilized  after  100  s,  which  is  not  acceptable.  Also,  the  current 
was  stabilized  at  a  value  0.1  A  below  the  desired  value  of  30  A. 
After  the  cell  voltage  decrease,  the  fuel  outlet  became  almost  com¬ 
pletely  starved  of  hydrogen  and  carbon  monoxide.  This  could  lead 
to  damage  of  the  tube.  Often  the  fuel  channel  outlet  stream  is 
combusted  in  a  tail  gas  combustor  to  provide  required  heat  for 
the  fuel  cell  system.  The  lean  outlet  condition  experienced  dur¬ 
ing  the  simulation  could  cause  the  flame  in  tail  gas  combustor  to 
extinguish. 

5.  Linear  identification  of  the  SOFC  stack 

The  physical  model  contains  the  coupled  effects  of  fuel  channel 
flow,  porous-media  transport,  heat  transfer,  reforming  chemistry, 
and  electrochemistry.  Although  already  simplified  compared  to 
even-more-complex  physical  models,  the  required  computational 
resources  are  still  too  great  for  direct  incorporation  into  an  MPC 
implementation  [14].  Although  the  physical  model  is  highly  non- 
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Fig.  9.  Total  current,  exhaust  H2  mole  fraction,  and  temperatures  with  respect  to 
time  in  going  from  a  low  power  demand  to  a  high  power  demand.  The  plots  capture 
tube  condition  changes  due  to  thermal  dynamics  which  are  on  the  order  of  hun¬ 
dreds  of  seconds.  The  results  shown  for  the  first  10  s  of  simulation  are  exactly  the 
same  as  those  illustrated  in  Fig.  8.  The  voltage  is  dropped  from  0.78  V  to  0.69  V  one 
second  after  the  start  of  simulation.  The  velocity  of  the  fuel  inlet  is  increased  from 
24.7  cms-1  to  51.4cms_1  over  half  a  second  starting  at  time  equal  to  1.2  s.  Over  the 
course  of  half  a  second,  the  air  flow  is  raised  from  33.5  mg  s-1  to  43.5  mg  s-1  starting 
80  s  into  the  simulation. 


linear,  it  is  possible  to  extract  low-order  locally  linear  models 
that  capture  the  dominant  dynamics  of  the  high-order  physi¬ 
cal  model  at  particular  OPs.  The  companion  paper  [1]  extends 
the  method  to  nonlinear  system  identification  over  a  range  of 
OPs. 

Model  reduction  is  accomplished  using  a  data-based  approach 
that  is  based  on  the  subspace  class  of  system  identification 
methods.  In  this  approach,  the  physical  model  takes  the  role 
of  an  experiment,  albeit  one  without  measurement  noise.  This 
approach,  which  is  an  alternative  to  the  direct  mathematical  lin¬ 
earization  of  the  physical  model,  has  several  advantages.  First, 
there  is  no  need  to  for  the  physical  model  to  be  represented 
as  ordinary  differential  equation  (ODEs)  in  standard  form  [i.e., 
y'  =  /(t,y)].  The  physical  model  used  here  is  expressed  as  DAEs, 
not  ODEs.  Second,  a  data-based  model  provides  the  best  lin¬ 
ear  approximation  of  the  nonlinear  system  as  measured  over 
the  amplitude  and  frequency  range  of  the  experimental  input, 
whereas  a  linearized  model  has  no  corresponding  approximation 
qualities. 

The  linear  identification  of  the  (SOFC)  stack  proceeds  as  fol¬ 
lows.  First,  a  nominal  OP  for  the  stack  is  selected,  (u,y),  where 
u  are  the  nominal  inputs  and  y  are  the  nominal  outputs  of  the 
stack.  A  small-signal  sequence  8u  is  designed  with  frequency  con¬ 
tent  that  matches  the  expected  system  bandwidth.  One  input  is 
perturbed  around  its  OP,  while  fixing  other  inputs  at  their  nomi¬ 
nal  values.  The  small-signal  response  8y  =  y  -  y  is  recorded.  This 
procedure  is  repeated  for  all  system  inputs.  From  this  data,  a 
single-input  single-output  (SISO)  reduced-order  model  of  each 
input-output  pair  is  identified.  These  SISO  models  are  combined 


into  a  single  multiple-input  multiple-output  (MIMO)  system,  and 
reduced  using  balanced  model-reduction  methods  [15],  resulting 
in  a  single  (MIMO)  reduced-order  linear  model  at  a  particular  OP. 
As  explained  briefly  in  the  following  section,  the  subspace  sys¬ 
tem  identification  is  used  to  capture  the  small-signal  model  at  a 
particular  OP. 

5.1.  Subspace  system  identification 

The  approach  to  fit  a  low-order  model  to  data  (here,  output 
from  the  physical  model)  is  based  upon  subspace-identification 
methods.  Although  this  process  is  well  known  and  documented 
(e.g.,  [16]),  it  is  mathematically  complex.  Only  a  brief  summary  is 
provided  here. 

An  input-output  time  sequence  ( 8uk,8yk\  (fc  =  l,...,N) 
obtained  and  recorded  by  exercising  the  physical  model  at  a  sam¬ 
pling  rate  Ts.  The  objective  is  to  find  a  linear  time-invariant  (LTI) 
system  in  state-space  form, 

/  Xfc+1  =  +  C8uk 

\ Syk  =  Cxk  +  DSuk  ’  11  ] 

that  can  closely  reproduce  the  data.  The  state-space  vector 
is  xkeRn  and  all  other  vectors/matrices  are  assumed  to  have 
compatible  dimensions.  The  LTI  model  (Eq.  (43))  expresses  the 
relationships  among xfc+1  ,xk,uk,  and yk  over  a  single  sampling  inter¬ 
val.  The  subspace  identification  proceeds  by  collecting  together 
larger  sampling  intervals  (windows)  of  the  input-output  data. 
The  window  length  s  is  a  free  parameter  that  should  be  chosen 
to  be  larger  than  the  expected  order  of  the  identified  system, 
n.  Considering  an  s-length  segment  of  data,  Eq.  (43)  can  be 
rewritten  as 
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Because  the  model  to  be  identified  is  time-invariant  (i.e., 
the  0,  r ,  C,  and  D  matrices  are  time-independent),  any 
shifted  window  of  the  data  can  be  expressed  as  Eq.  (44).  The 
shifted-window  input-output  data  can  be  expressed  in  matrix 
form  as 


Ys  =  rsX  +  HsUs, 
where 


(45) 
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Fig.  10.  Small-signal  perturbation  with  fuel  flow  rate  changing  with  sampling  time 
of  Ts  =0.125  (s). 


The  matrix  rs  is  called  the  extended  observability  matrix  where 
s  >  n  guaranties  rs  be  full  column  rank  (observability  condition). 
Particular  system  parameters  are  contained  in  Hs  and  rs,  and  can 
be  extracted  when  these  terms  are  known.  Following  the  subspace 
system-identification  approach,  estimates  of  0,  r,  C,  and  D  matri¬ 
ces  can  be  evaluated  [17].  These  steps  are  briefly  explained  here. 
The  matrices  Ys  and  Us  (Eq.  (45))  are  known  from  the  “measure¬ 
ments”  obtained  by  exercising  the  physical  model  with  small-signal 
inputs  Su  and  recording  the  outputs  by.  The  HSUS  term  can  be 
removed  by  multiplying  Eq.  (45)  from  right  by 

=  /  -  uj(usuj  )_1  us,  (50) 


where  Y1ut±  is  the  orthogonal  complement  of  Us.  As  a  result,  Eq. 
(45)  becomes 

YsnUT±  =  rsxnUT±.  (5i) 


The  singular  value  decomposition  (SVD)  of  Ysnffr±  can  be  rep- 

us 

resented  as 

(52) 

where  5  and  V  are  orthonormal  matrices  and  £  is  a  diagonal  matrix 
containing  the  singular  values.  These  matrices  can  be  further  par¬ 
titioned  in  two  parts  as 

SEI/t  =  [S,S2] 

where  £i  contains  the  first  n  dominant  singular  values.  Under 
appropriate  condition  (i.e.,  persistently  exciting  input),  the  range 
space  of  Si  is  the  same  as  rs.  Thus,  an  estimate  for  rs  is  taken  to 
be  the  n  first  columns  of  S  (Fs  is  defined  only  within  a  coordinate 
transformation  ofxk).  With  rs  known,  the  matrix  C (Eq.  (43))  is  the 
first  p  rows  of  rs ,  where  p  is  the  number  of  system  outputs.  The 
matrix  0  can  be  determined  through  the  shift  pattern  of  rs.  With 
0  and  C  in  hand,  r  and  D,  which  are  linear  in  Eq.  (45),  can  be  found 
via  solving  a  linear  least  squares  problem. 

5.2.  Illustration  of  linear  system  identification 

This  section  illustrates  the  results  of  linear  identification  for  the 
(SOFC)  stack  at  a  particular  OP.  The  (SOFC)  stack  variables  are  par¬ 
titioned  into  three  inputs  and  four  outputs.  The  input  variables  are 
cell  voltage,  fuel  mass  flow  rate,  and  air  mass  flow  rate.  The  output 
variables  are  cell  current,  hydrogen  concentration  in  exhaust,  aver¬ 
age  MEA  temperature,  and  cathode-exhaust  air  temperature.  As 
described  earlier,  the  physical  model  represents  a  widely  disparate 
range  of  characteristic  time  scales.  Electrical  response  is  essentially 
instantaneous  (i.e.,  because  electrochemical  double-layer  charging 
is  neglected,  a  change  in  operating  voltage  causes  an  instantaneous 
change  in  current).  Characteristic  response  time  for  fluid  flow  and 
diffusion  is  of  the  order  of  one  second  or  less.  The  characteristic 
times  for  thermal  response  is  much  longer,  on  the  order  of  min¬ 
utes.  For  the  purposes  of  illustrating  the  system  identification,  air 
flow  rate  and  cell  temperatures  are  considered  constant.  Thus,  the 
illustration  here  focuses  upon  identifying  cell  current  and  hydrogen 
concentration  from  cell  voltage  and  fuel  flow  rate.  In  the  compan¬ 
ion  paper  [1  ],  a  system  is  identified  for  the  full  range  of  times  scales 
and  an  MPC  controller  is  designed. 

Assume  a  nominal  OP  for  a  stack  that  is  operating  at  a  volt¬ 
age  of  Eceu  =  0.74  V,  fuel  mass  flow  rate  of  5.2  mgs-1,  and  air  mass 
flow  rate  of  33.5  mg  s-1 .  Under  these  conditions,  the  physical  model 
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Fig.  11.  Comparison  between  low-order  and  high-order  models  for  validation  simulation. 
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(single  tube)  predicts  a  net  current  of  24.96  A  with  10.24%  H2 
mole  fraction  in  fuel  outlet.  The  stack  is  excited  by  a  perturbation 
sequence  that  is  suitable  for  process  model  identification.  In  this 
example,  a  pseudo  random  binary  sequence  (PRBS)  with  zero  mean 
value  is  used  as  the  perturbation  sequence  [18,19].  A  sequence  of 
500  samples  is  used,  with  each  input  perturbed  separately.  Fig.  10 
shows  a  segment  of  a  simulation  at  the  selected  OP.  In  this  per¬ 
turbation,  the  fuel  flow  rate  is  perturbed  with  a  sampling  time  of 
Ts  =  0.125  s,  while  the  other  inputs  are  held  constant  at  their  nom¬ 
inal  operating  values.  Fig.  1 0a  shows  the  PRBS  that  is  applied  as  the 
fuel  flow  variation.  Fig.  10b  shows  the  corresponding  change  in  the 
H2  mole  fraction  in  fuel  exhaust.  The  subspace  identification  at  this 
OP  results  in  a  12th-order  small-signal  model. 

To  validate  the  small-signal  model,  an  input  sequence  is 
designed  in  which  all  the  inputs  vary  simultaneously.  The  identified 
small-signal  model  should  be  able  to  predict  the  output  of  this  sim¬ 
ulation.  Fig.  11a  and  b  shows  the  input  variations.  Fig.  11c  shows  the 
predicted  cell  current,  comparing  the  high-order  physical  model 
and  the  identified  low-order  linear  model.  In  this  example,  the  two 
curves  are  nearly  indistinguishable.  Fig.  lid  shows  another  vali¬ 
dation  result  comparing  the  hydrogen  mole  fraction  in  the  anode 
exhaust.  Again,  the  comparison  is  excellent. 

6.  Conclusion 

A  physically  based,  transient  model  for  a  tubular  anode- 
supported  (SOFC)  is  developed  as  the  basis  for  implementing  an 
MPC  controller.  Considering  the  coupled  interactions  of  fuel  flow, 
porous-media  transport,  heat  transfer,  reforming  chemistry,  and 
electrochemistry,  the  model  is  implemented  with  error-controlled 
DAE  software  that  enables  the  accurate  prediction  of  accurate  tran¬ 
sient  responses.  Accounting  for  widely  disparate  time  scales,  the 
model  resolves  spatial  and  temporal  profiles  of  composition,  tem¬ 
perature,  current  density,  and  velocity  throughout  the  cell.  Even 
with  approximations  to  reduce  the  complexity  of  the  physical 
model,  it  is  still  to  large  to  be  incorporated  directly  into  an  MPC 


implementation.  Thus,  linear  model  reduction  is  used  to  reduce  the 
high-order  physical  model  to  a  low-order,  locally  linear  model.  An 
example  is  used  to  illustrate  the  validity  of  the  reduced  model  at  a 
particular  OP.  As  discussed  in  a  companion  paper  [1  ],  the  reduced- 
order  models  form  the  basis  for  implementing  an  MPC  controller. 
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